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ABSTRACT 

A Mathematical Model of a Copper Oxide/Copper “Vaporvolt’ Cell. (August 1992) 
Makoto Kawanami, B.Eng., Kyoto University; 

M.Eng., Kyoto University 
Chair of Advisory Committee: Dr. Ralph E. White 

A new type of battery named “Vaporvolt”* cell is in the early stage of its 
development. A mathematical model of a CuO/Cu “Vaporvolt” cell is presented 
that can be used to predict the potential and the transport behavior of the cell 
during discharge. A sensitivity analysis of the various transport and electrokinetic 
parameters indicates which parameters have the most influence on the predicted 
energy and power density of the “Vaporvolt” cell. This information can be used to 
decide which parameters should be optimized or determined more accurately through 
further modeling or experimental studies. The optimal thicknesses of electrodes and 
separator, the concentration of the electrolyte, and the current density are determined 
by maximizing the power density. These parameter sensitivities and optimal design 
parameter values will help in the development of a better CuO/Cu Vaporvolt cell. 


* Product name of Energy Innovations, Inc. 


IV 


TABLE OF CONTENTS 


ABSTRACT 

LIST OF TABLES 

LIST OF FIGURES 

INTRODUCTION 

MODEL DEVELOPMENT 

Electrochemical Reactions at Electrodes 

Assumptions 

Governing Equations 

Cell Voltage 

RESULTS AND DISCUSSION 

The Effects of the Separator Thickness 
The Effects of Current Density .... 

Sensitivity Analysis 

Optimization 

SUMMARY 

NOMENCLATURE 


Page 
. iii 
. vi 
. vii 
. . 1 
. . 4 
. . 4 
. . 5 
. . 6 
. 11 
. 16 
. 16 
. 17 
. . 17 
. . 29 
. . 32 
. . 33 


REFERENCES 


35 



V 


APPENDIX 37 

Computer Program Data File 37 

Computer Main Program File 3 ^ 

VITA 54 



VI 


LIST OF TABLES 

Page 

Table I. Base case parameter values for the CuO/Cu “Vaporvolt cell 

simulation 

Table II. Comparison of the performances of the “Vaporvolt” cells 

with different separator thicknesses 16 

Table III. The result of the optimization 29 

Table IV. Comparison of the “Vaporvolt” cell performances with 

lithium batteries 50 



vn 


LIST OF FIGURES 

Page 

Figure 1. Schematic of a “Vaporvolt” cell 7 

Figure 2. Cell voltage vs. capacity delivered 18 

Figure 3. The effects of current density on capacity of “Vaporvolt” 

cells 

Figure 4. The effects of current density on energy density of 

“Vaporvolt” cells 20 

Figure 5. The effects of current density on power density of 

“Vaporvolt” cells 21 

Figure 6. Sensitivity of the model predictions to the changes in the 

current density 23 

Figure 7. Sensitivity of the model predictions to the changes in the 

concentration of KOH 24 

Figure 8. Sensitivity of the model predictions to the changes in the 

thickness of the separator. 25 

Figure 9. Sensitivity of the model predictions to the changes in the 

thickness of the positive electrode 26 

Figure 10. Sensitivity of the model predictions to the changes in the 

thickness of the negative electrode 27 


The effects of current density on power density of the cell 
when thicknesses of the electrodes and separator are 
optimized 



1 


INTRODUCTION 

Development of batteries with higher energy and power densities has been desired 
for many years. The traditional approach has been to employ more energetic reactions, 
which results in higher unit cell voltage, to achieve higher power and energy densities 
(1-13). Although stored energy and power densities have both increased by this 
approach (14), the specific stored energy and power densities are still relatively low. 
This is because of the additional weight required to contain these highly reactive 
materials and to maintain the mechanical strength of the battery. In addition, the 
energy densities of secondary batteries are much lower than that of primary batteries 
(15). Furthermore, there is a safety concern associated with using highly reactive 
chemicals such as lithium and thionyl chloride (16-18). Another approach to achieve 
higher specific energy and power densities is to minimize unnecessary weight and 
volume in existing systems. 

“Vaporvolt” cells, which are currently being evaluated by Energy Innovations, Inc. 
for NASA as possible replacements for high power density primary lithium battenes, 
use this approach. Weight and volume are minimized by using thin laminations of 
metals and their oxides as the active materials of a battery (19). Reaction rates at 
the electrode surfaces of these batteries are low, resulting in minimal stress on the 


This document follows the style of the Journal of the Electrochemical Society. 
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electrodes. High currents are obtained by using the high area per volume character- 
istic of thin laminations. Desired voltage can be obtained by bipolar construction. 
Electrochemical conversion is accomplished by displacing oxygen from one metal to 
the other. There are no side reactions that complicate the process and make unwanted 
by-products; and, in most cases, the energy conversion process is highly reversible. 

In designing high performance “Vaporvolt” cells, there are various attributes that 
can significantly influence the system. Such attributes might be the thicknesses of the 
electrodes and the separator, KOH concentration, and the current density. A sensitivity 
analysis is performed on various parameters to determine which parameters are most 
influential in increasing or decreasing the energy and power densities. This information 
can indicate the direction one should take in order to design better Vaporvolt cells. 
The results of the sensitivity analysis can also suggest which parameters should be 
obtained with more accuracy through further modeling studies or experimentation. 

To achieve high performance in the “Vaporvolt” cell, various design parameters 
can be optimized so that the cell will deliver the maximum attainable power density. 
The important parameters in the “Vaporvolt” cell are the thicknesses of the positive 
electrode ( L pos ), negative electrode (Lmg) and separator (L s ), the concentration of KOH 
( c KOh)i and current density (i). The CuO/Cu system was selected in this investigation 
because the CuO/Cu “Vaporvolt” cells are the only successfully manufactured batteries 
at this time. By using a mathematical model of the CuO/Cu “Vaporvolt” cell, these 
parameters are investigated in order to determine if an optimal value exists for each 
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parameter. The sensitivity of the model predictions to various parameters will then be 
examined followed by the determination of the optimal design parameters to maximize 
the power density of the “Vaporvolt” cell. 
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MODEL DEVELOPMENT 

Electrochemical Reactions at Electrodes 

The reactions occurring in the CuO/Cu “Vaporvolt” cell during discharge are the 
reduction of CuO at the cathode (positive electrode) 

Positive Electrode : CuO + — H 2 O + e — ► — CU 2 O + OH [1] 

Upos = 0.669 — 0.0591 pH V vs. NHE [2] 

and the oxidation of copper at the anode (negative electrode) 

Negative Electrode : Cu + OH — » — CU 2 O + — H 2 O + e [3] 

Uneg = 0.471 - 0.0591pH V vs. NHE [4] 

where U pos and U neg are the open circuit potentials of the positive and negative 
electrodes, respectively. The electrochemical data in equations [2] and [4] are taken 
from reference (20), at 25 °C. 

The overall reaction is the production of CU 2 O from CuO and Cu. 


Total : 


CuO + Cu -+ Cu 2 0 


[5] 


Since the solubility of Cu 2 0 in basic solution is very low, solid Cu 2 0 is formed 
during discharge. The open circuit voltage of a CuO/Cu ‘ Vaporvolt cell, U ce ii is 
obtained by subtracting equation [4] from equation [2]. 
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U ce u = 0.198 V vs. NHE 


[ 6 ] 


This open circuit voltage of CuO/Cu couple is not so high. It is easily predicted that 
if we employ a system with higher open circuit voltage, we can obtain higher energy 
and power density. However, CuO/Cu couple is the only successfully manufactured 
“Vaporvolt” cells at this time. 

The objective of the present investigation is to develop a mathematical model 
that could be used to estimate the optimum values of transport and electrokinetic 
parameters, and to predict the maximum attainable energy and power densities of a 
CuO/Cu “Vaporvolt” cell. 

Assumptions 

The assumptions are as follows. 

1. One dimensional model can describe the behavior of the cell. 

2. Dilute solution theory (21), applies. Negligible interactions between the solutes 


are assumed. 
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3. The Nemst-Einstein equation, D, = mRT, which is implicit in the dilute solution 
theory, applies. 

4. The current density takes the form of the Butler- Volmer equation, which expresses 
the exponential dependence of the current on the overpotential. 

5. The physical, transport, and electrokinetic parameters are constant throughout the 
solution. 

6. The cell is isothermal. 

7. No homogeneous chemical reactions occur in the electrolyte. 

8. No electrolyte movement. 

9. The electrochemical reactions occur at the surface of the electrodes. 


Governing Equations 

A schematic of a CuO/Cu “Vaporvolt” cell is shown in Fig. 1. During 
discharge of a CuO/Cu “Vaporvolt” cell, hydroxide ions are released from the positive 
electrode (see equation [1]), and consumed at the negative electrode (see equation [3]). 
Therefore, hydroxide ions must be transported in the electrolyte in the separator from 
the positive electrode to the negative electrode. Mass transport in this system is due 
to migration in an electric field and diffusion in a concentration gradient. Therefore, 
the flux expression for each species i can be written as 


N, = — z,u,Fc,V$ — Z),Vcj 


[7] 
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where the ionic mobility, u,, is assumed to be related to the diffusion coefficient 
Di by the Nemst-Einstein equation 


u, = 


D \ 
RT 


[ 8 ] 


The first and the second terms of right side of equation [7] represent the mass transport 
in solution due to migration in an electric field and diffusion in a concentration 
gradient, respectively. The differential material balance equation for each species 
without homogeneous chemical reaction in the electrolyte is as follows 


dc t 

dt 


= -V-N, 


[9] 


By combining equations [7] and [9], we can obtain following equation 


^ = z t u,F(c,V 2 4> + Vc t V$) + AV 2 c, 


[ 10 ] 


Since the net charge in the bulk electrolyte is zero, the equation of electroneutrality 



I 


[ 11 ] 
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must hold in the electrolyte. 

The flow of charge is related to the current density i in the solution: 


i = F^ziNi 
i 


[ 12 ] 


Equations [7] and [12] can be combined and rearranged to get 


where 


V$> = — - — — ( Y ziDiVci 

K K \ ' 


* = h V zlDiCi 


[13] 


[14] 


Equations [10] and [13] for one spatial coordinate reduces to equations [15] and 
[16], respectively 


dc t 

Hi 


z{u,F 



d*d$\ 

^ dx dx ) 


+ Di 


d 2 c, 

dx 2 


[15] 


d<& -i F ( dci \ 


[16] 


There are three unknowns ck+ , cqh- and 4>, so we need three governing equations. 
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Equation [11] is used for c K +, equation [15] is used for c 0 h~> and equation [16] is 
used for 

The initial and the boundary conditions for these dependent variables are as 
follows. 

At f=0, the concentration of the species are their initial values 

For all x at t = 0 c t (t = 0, x) = c° [17] 


since c° must obey Y^ z i c ° ~ 0- 

t 

At x=0 (positive electrode surface), the change in the concentration of OH' is related 
to the current density since the number of OH' ions released at the positive electrode 
surface is proportional to the number of the electrons received at the electrode, that 
is, current density (see equation [1]). 


For all t, t > 0 


dean- . _ j 

dx x ~° ti\FDoh- 


[18] 


where n\ is the number of the electrons transferred in equation [1]. 

At x=L s (negative electrode surface), the change in the concentration of OH' is also 
related to the current density (see equation [3]). 

dcoH- | 

\x=L, = 


For all t, t > 0 


dx 


n^FDoH- 


[19] 
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where nj is the number of the electrons transferred in equation [3]. 

Since we cannot measure the absolute potential, we need to fix a potential at x = 0: 

For all t at x = 0 <!>(<, x = 0) = 0 [20] 


Cell Voltage 

The relation between the current density, /, and the overpotential, r], is expressed 
by the Butler- Volmer electrochemical rate expression 



where A* and A are the active surface area of the electrode and the geometric surface 
area of the electrode, respectively. Active surface area of the electrodes decreases 
with the discharge of the cell since Cu 2 0, that has very low electric conductivity, is 
formed on the electrodes during discharge. The active surface area of each electrode 
is assumed to change according to the following equation 


A* 



[ 22 ] 


where m is the amount of active material (CuO for positive electrode, Cu for negative 
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electrode) in an electrode at time t, and m° is the initial amount of active material 
in the electrode. The initial amount of active material, m°, is simply expressed by 
the following equation 


m° = Ljdj 


[23] 


where Lj is the thickness of the electrode and dj is the density of the electrode. 
For a constant current density discharge, the amount of active material, m, changes 
with time according to following equation 


m = m° 


Mj i t 
rijF 


[24] 


where Mj is a molecular weight of an electrode j. 

By combining equations [21] through [24], we can obtain following relation between 
the current density i and the overpotential t] 



Mj i t 

LjdjUjF 



[25] 


Note that CuO, Cu 2 0 and Cu exist as solids because of their low solubility in a 
basic solution, so the values of ar *d Cu 2 0 and Cu in equation 
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[25] remain constant. 

In addition, the anodic and cathodic transfer coefficients are assumed to sum to the 
number of electrons transferred (i.e., a a + a c — rij), and the reaction order coefficients 
are assumed to be simply related to s ,• 


Pi = S{ qi = 0 if s t > 0 

Pi = 0 q, = —Si if s,< 0 


[26] 


And the electrode potentials, E pos and E neg , are expressed by 


Epos — Epos “t” — fi) "T Vpos 


[27] 


Eneg — Eneg j ® — L $ ) "t" Tjneg 


[28] 


The potential of a “Vaporvolt” cell, E ce p, is the difference between the potentials 
of the positive electrode and the negative electrode. By subtracting equation [28] 
from equation [27], 


E ce H — E pos E, 


-'neg 


— Ecell ~ i % — L s ) + T] p os Ineg 


[29] 


U C ell — Epos E n eg 


[ 30 ] 
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where U ce u is the open circuit potential of the “Vaporvolt” cell. 

Note that an overpotential takes a negative value at the positive electrode, and 
takes a positive value at the negative electrode, and $(<,x = L s ) is always positive 
during discharge because negatively charged particle (OH - ) must transport from the 
positive electrode ( x = 0) to the negative electrode (x = L s ). Therefore, the potential 
of a “Vaporvolt” cell, £«.//, is always smaller than the open circuit potential of the 
“Vaporvolt” cell, U ce u- 

Equations [7] - [16] represent the basic governing equations necessary to describe 
the behavior in the electrolyte of a “Vaporvolt” cell. It should be noted that these 
equations are solved numerically by setting the cell current density, i, and calculating 
the potential of the “Vaporvolt” cell, E ce u, or the power density, P (= iE ce u). B AND(J), 
a finite difference numerical technique developed by Newman (21), was used to solve 
the equations. The computer program used for this investigation is listed in Appendix. 
The model parameters associated with these equations are shown in Table I. 
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Table I. Base case parameter values for the CuO/Cu “Vaporvolt' 


Structural parameters: 

Thickness of positive electrode (CuO) ( L pos ) 

Thickness of negative electrode (Cu) ( L neg ) 

Separator thickness ( L s ) 

Open circuit potential ( U ce u ) 

Current density (i) 

Anodic and cathodic exchange current densities (i°) 
Anodic transfer coefficient for Cu oxidation (a a ,Cu) 
Cathodic transfer coefficient for CuO reduction (a Cj cuo) 
Electrolyte concentration {etc oh) 

Initial K + concentration (c°^ + ) 

Initial OH - concentration {c° OH ~) 

Reference K + concentration (c r K+ ) 

Reference OH - concentration ( c T OH _) 

Temperature (T) 


cell simulation. 


0.00254 cm 
0.00254 cm 
0.00254 cm 
198 mV 
0.001 A/cm 2 
0.001 A/cm 2 
0.5 
0.5 

8.3878 M 
8.3878 M 
8.3878 M 
1.0000M 
1.0000M 
298 K 
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RESULTS AND DISCUSSION 

The Effects of the Separator Thickness 

It is assumed that the electrochemical reaction occurs only at the surface of the 
active material, however; actual reaction layer moves into the inside of the active 
material as the cell discharges. Therefore, the thickness of the mass transport region 
of the cell increases as the cell discharges. To check the validity of this assumption, 
the cell performances with different separator thicknesses are predicted and compared. 
First, the cell performance is predicted with the parameter values listed in Table I, and 
then, the thickness of the separator is doubled, and the cell performance is predicted 
again. The result is shown in Table II. The predicted performances of “Vaporvolt” 
cells with different separator thicknesses are almost the same. This result indicates 
that the movement of the reaction layer in the cell has little effect on the performance 
of the cell. Therefore, the assumption that the electrochemical reaction occurs only at 
the surface of the active material is valid in this model. 

Table II. Comparison of the performances of the “Vaporvolt” 
cells with different separator thicknesses. 


Separator thickness 

0.00254 cm 

0.00508 cm 

Cell capacity 

17.82711 C/cm 2 

17.82711 C/cm 2 

Energy density 

1.88573 Ws/cm 2 

1.88569 Ws/cm 2 

Power density 

0.10637 mW/cm 2 

0.10636 mW/cm 2 


17 


The Effects of Current Density 

Fig. 2 shows the cell voltage change during the discharge with the parameter 
values listed in Table I, at various current densities. The cell voltage drop due to 
the increases in the overpotentials at the positive and the negative electrodes as we 
increase the current density during discharge. Therefore, the capacity and the energy 
density of the cell decrease with increasing current density of the cell during discharge, 
as shown in Figs. 3 and 4. However, the power density of the “Vaporvolt” cell 
is not a monotonic function of the current density. As shown in Fig. 5, the power 
density of the cell first increases and then decreases as the current density increases. 
So there exists a maximum power density value at a current density. 


Sensitivity Analysis 

In order to determine the relative importance of the design feature of the cell 
(thickness of the electrode, e.g.) and the transport and kinetic parameters on the 
“Vaporvolt” cell’s performance, a sensitivity analysis was done. The sensitivity 
analysis can indicate which parameters have the largest influence on the predicted 
performance of the “Vaporvolt” cell. If a small perturbation in a parameter does not 
significantly change the predicted capacity, energy density or the power density of 
the “Vaporvolt” cell, then that parameter could assume a large range of values, all of 
which will give the same performance. The sensitivity coefficient can be defined as the 
dimensionless change in the predicted capacity, energy density, and the power density 
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Figure 2. Cell voltage vs. capacity delivered. 
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Figure 3. The effects of current density on capacity of “Vaporvolt” cells. 
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Figure 5. The effects of current density on power density of “Vaporvolt” cells 
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of the cell for a small dimensionless perturbation in a parameter^', while holding all 
other parameters constant (22): 


Sensitivity = 


dG/G 

dOj/Oj 


[31] 


where 6 } is the parameter and G is the capacity, the energy density, or the power 
density of the cell. 

Hence, large sensitivity coefficients mean that the parameter of interest signifi- 
cantly influences the cell performances. If the sensitivity coefficients for a parameter 
is large, the parameters should be obtained with more accuracy through further mod- 
eling or experimental studies. That is, if the value for a parameter is not accurately 
known and the parameter has a large sensitivity coefficient, then that parameter value 
should be ascertained more accurately to gain confidence in the model predictions. 

All sensitivity coefficients calculated for this work were accomplished by increas- 
ing the parameter of interest by 5% over the base case value (shown in Table I) and 
calculating the resulting change in the capacity, energy density, and power density of 
the cell from the base case values. The sensitivity coefficients for various parameters 
are shown in Figs. 6 through 10. As shown in Fig. 6, the model predictions are most 
sensitive to the current density. 

The model predictions show little sensitivity to small perturbations in the KOH 
concentration as shown in Fig. 7. This is because the IR drop due to the resistance 



Figure 6. Sensitivity of the model predictions to the changes in the current density 
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Figure 7. Sensitivity of the model predictions to the 
changes in the concentration of KOH. 





I 


-2 


10 _1 


1 


Separator thickness 


/ cm 


Figure 8. Sensitivity of the model predictions to the 
changes in the thickness of the separator. 
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Figure 9. Sensitivity of the model predictions to the changes 
in the thickness of the positive electrode. 
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Figure 10. Sensitivity of the model predictions to the changes in 
the thickness of the negative electrode. 
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of the electrolyte is negligible since the thickness of the separator is very small. The 
concentration of KOH in this region has very small effect on the capacity, the energy 
density, and the power density of the cell. 

The effects of the “Vaporvolt” cell separator thickness on the predicted perfor- 
mance are shown in Fig. 8. The sensitivity coefficients for the capacity, the energy 
density, and the power density of the cell decrease similarly and monotonically with 
the increase of the separator thickness. Therefore, the thickness of the separator should 
be as thin as possible, however; to avoid shortening of the cell, the thickness of the 
separator should be about 0.001 cm. 

Figs. 9 and 10 show the effects of “Vaporvolt” cell positive and negative electrodes 
thickness on the predicted performances. The capacity and the energy density of the 
“Vaporvolt” cell increase as we increase the electrode thicknesses up to 0.002 cm since 
the sensitivity coefficients for the capacity and the energy density are positive in this 
region, but the capacity and the energy density decrease as we increase the electrode 
thicknesses when the thicknesses are above 0.002 cm. Therefore, the capacity and 
energy density become their maximum when the electrode thicknesses are about 0.002 
cm. The sensitivity coefficient for the power density of the cell is always negative, as 
shown in Figs. 9 and 10, so a small electrode thickness is favorable to get a larger 
power density. However; employing electrodes that are less than 0.001 cm could 
results in shortening of the cell. A thickness of 0.001 cm is favorable for the positive 
and negative electrodes of a “Vaporvolt” cell to obtain maximum power density. 
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Optimization 

The maximum power density can be obtained when we employ 0.001 cm as the 
thicknesses of the separator, the positive electrode, and the negative electrode, and 
current density 2.5 mA/cm*. As shown in Table 111, the maximum power density 
obtained is 45.1 W/l, and the capacity is 1768 kC/1 and the energy density is 26.5 
Wh/1 at these parameter values. These values are lower than that of the htgh 
energy density primary cells, such as Li/SOCl 2 cells, which highest power density 
is approximately 80 W/l and its energy density is about 150 Wh/1 when the highest 
power density is obtained (14), however, the maximum attainable power densrty of a 
CuO/Cu “Vaporvolt” cell is much higher than that of Li/CF, or Lr/MnCh cells, which 
highest power density is approximately 15 W/l and the energy density is about 50 
Wh/1 when the highest power density is obtained (14). This is summarized in Table 
IV. If stable electrodes with thickness less than 0.001 cm could be developed in the 
future, much larger power density could be obtained. For example, if we could make 

Table III. The result of the optimization. 
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electrodes and separators with thicknesses of 0.0005 cm, the maximum power density 
is 90.3 W/l, as shown in Fig. 11; and the capacity and the energy density wxll be 

1768 kC/1 and 26.5 Wh/1, respectively. 


Table IV. Comparison of the “Vaporvolt” cell performances with lithium batteries. 


CuO/Cu 

"Vaporvolt" 


Li/SOCh Li/Mn0 2 


Power density, W/l 
Energy density, Wh/1 
Unit cell open circuit voltage , V 


45.1 

80* 

15 

26.5 

150* 

50* 

0.198 

3.67* 

3.5* 


* These values are taken from ref. [14]. 



Optimum value 

Figure 11. The effects of current density on power density of the cell 
when thicknesses of the electrodes and separator are optimized. 
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SUMMARY 

A mathematical model of a CuO/Cu “Vaporvolt” cell is presented to predict the 
potential and the transport behavior of the cell during discharge. A sensitivity analysis 
of a CuO/Cu “Vaporvolt” cell model indicates that several parameters can significantly 
influence the performance of the system. In particular, current density has been found 
to influence significantly the performance. The effect of various design parameters 
have been investigated to determine if optimal values exist for the parameters. The 
model has shown that the thicknesses of the electrodes and separator, and current 
density can be optimized to give the maximum attainable power density. The model 
has also shown that the power density can be significantly improved if we could make 
thinner electrodes and separators could be made. 
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NOMENCLATURE 


Roman Symbols 

/C Active surface area of the electrode, cm 2 

A Surface area of the electrode, cm 2 

d Concentration of species i, mol/cm 3 

Cf° Initial concentration of species mol/cm 

3 

c{ Reference concentration of species /, mol/cm 

Di Diffusion coefficient for species /, cm 2 /s 

dj Density of electrode j, g/cm 3 

E Electrode potential, V 

F Faraday’s constant, 96487 C/mol 

G Capacity, Energy density, or Power density 


i 

i 

i° 

L 

Mj 

m 

m° 


N,- 


P 

Pi 

Qi 

R 


Ri 


Si 

T 


Current density vector, A/cm 2 
Current density, A/cm 2 
Exchange current density, A/cm 2 
Thickness, cm 

Molecular weight of electrode j, g/mole 
Amount of active material at time t, g/cm 2 
Initial amount of active material, g/cm 2 
Number of electrons transferred in reaction j 

Flux vector of species /, mol/(cm 2 -s) 

Power density, W/l 

Anodic reaction order for species i 

Cathodic reaction order for species i 

Gas constant, 8.3143 J/(mol-K) or 82.057 cm 3 atm/(mol-K) 

Production rate of species mol/(cm 3 -s) 

Stoichiometric coefficient of species / 

Temperature, K 
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Open circuit potential, V 
Mobility of species /, mol-cm 2 /(J-s) 
Electrolyte velocity vector, cm/s 
Charge number of species i 


Greek Symbols 


Anodic transfer coefficient 
Cathodic transfer coefficient 
conductivity of a solution, S/cm 
Overpotential, V 
Solution phase potential, V 
Parameter j 


Superscripts and Subscripts 

a Anode 


Cathode 
Species i 

Reaction j or electrode j 
Negative electrode 
Positive electrode 
Reference condition 
Separator layer 
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APPENDIX 


Computer Program Data File 


ishow 

0 

ishow=l for monitor output of 
iterations 

jmax 

100 

Max program iterations 

n j 

101 

Number of meshpoints 

nt 

400 

Number of timesteps 

ntpr 

40 

Print every ntpr timesteps 

nc 

6 

Number of current steps 

area 

4 . 00000000d+00 

Area of the separator, sq.cm 

cio 

8 . 38780000d-03 

Initial KOH concentration, mol/cc 

ecur (1) 

1 . 00000000d-03 

Exchange current density 
(positive electrode), A/cm2 

ecur (2) 

1 . 00000000d-03 

Exchange current density 
(negative electrode). A/ cm2 

aa (1) 

0 . 50000000d+00 

Anodic transfer coefficient (pos) 

ac (1 ) 

0 . 50000000d+00 

Cathodic transfer coefficient 
(pos) 

aa (2 ) 

0 . 50000000d+00 

Anodic transfer coefficient (neg) 

ac (2 ) 

0 . 50000000d+00 

Cathodic transfer coefficient 
(neg) 

ocv 

198 . 000000d-03 

Open circuit voltage, V 

delta 

2 . 54000000d-03 

Separator thickness, cm 

dpos 

2 . 54000000d-03 

Positive electrode thickness, cm 

dneg 

2 . 54000000d-03 

Negative electrode thickness, cm 

denp 

6 . 40000000d+00 

Positive electrode density, g/cm3 

denn 

8. 92000000d+00 

Negative electrode density, g/cm3 

denpd 

6 . 00000000d+00 

Positive electrode density 
after discharge, g/cm3 

dennd 

6 . 00000000d+00 

Negative electrode density 
after discharge, g/cm3 

atomp 

79 . 5454000d+00 

Atomic weight of positive 
active material, g/mole 

atomn 

63 . 5460000d+00 

Atomic weight of negative 
active material, g/mole 

atompd 

71 . 5430000d+00 

Atomic weight of p.a.m. 
after discharge, g/mole 

atomnd 

71 . 5430000d+00 

Atomic weight of n.a.m. 
after discharge, g/mole 

epos 

1 . 00000000d+00 

Number of electron transfer 



eneg 

1 . 00000000d+00 

diff (1) 

2 . 19000000d-05 

diff (2) 

2 . 19000000d-05 

temp 

298 . 000000d+00 

toler 

1 . 00000000d-05 

zi (1) 

-1 . 00000000d+00 

zi (2) 

1 . 00000000d+00 

cur (1) 

1 . 00000000d-03 

cur (2) 

2 . 00000000d-03 

cur (3) 

3 . 00000000d-03 

cur (4) 

4 . 00000000d-03 

cur (5) 

5 . 00000000d-03 

cur (6) 

6 . 00000000d-03 


(positive electrode) 

Number of electron transfer 
(negative electrode) 

Anion diff. coeff./ sq.cm/sec 
Cation diff. coeff., sq.cm/sec 
Temperature, K 
Relative error tolerance 
Anion (0H-) charge number 
Cation (K+) charge number 
Current density (1), A/cm2 
Current density (2), A/cm2 
Current density (3), A/cm2 
Current density (4), A/cm2 
Current density (5), A/cm2 
Current density (5), A/ cm2 
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Computer Main Program File 


c model2.f 

c last edition May. 10, 1992 

implicit real*8 (a-h,o-z) 

dimens ilpqon z ( 7 , 500) , a (7 , 7 ) ,b(7, 7) ,c(7,500) ,d(7, 15) , 

0 g(7),x{7,7),y(7,7), zk(500) , zi (7) , eta (2, 500) , 

& pot (500) , volt (500) , actp (500) , act n (500) , ecur (2) , 

& aa (2) , ac (2) , aarea (2) ,dif f (7) , cur (7) 

common A, B, C, D, G, X, Y, N, NJ 
OPEN (1, FILE='model2 . dat' , STATUS=' old' ) 

OPEN (2, FILE =, model2 . out' , STATUS=' unknown' ) 

OPEN (3, FILE=' model2. outt' ,STATUS=' unknown' ) 


c Model inputs. 

read (1, 9911) ishow, jmax, n j, nt, ntpr, nc 

9911 format (lOx, il5) 

read (1, 9912) area, cio, ecur (1) , ecur (2) , 

& aa (1) , ac (1) , aa (2) , ac (2 ) , ocv, delta, 

& dpos,dneg, 

& denp, denn, denpd, dennd, atomp, atomn, atompd, 

& atomnd, epos, eneg, dif f (1) , dif f (2) 

read(l, 9912) temp, toler, zi (1) , zi (2) 
read (1, 9912 ) (cur ( jc ) , jc=l, nc) 

9912 format (lOx, dl5 . 7) 
dposd=dpos*denp*atomp/ atompd/ denpd 

dnegd=dneg*denp*atomn/atomnd/dennd 

c Start current loop 

do 501 jc=l , nc 
write (6,721) jc,nc 

721 format (//, 10X, ' Current loop jc=',i3,' nc=',i3,/) 

do 720 i=2, 3 

write (i, 9913) ishow, jmax, n j, nt, ntpr 

write (i, 9914) cur ( jc) , area, cio, ecur (1) , ecur (2) , 

& aa ( 1 ) , ac ( 1 ) , aa (2) , ac (2) , ocv, delta, 

& dpos, dneg, dposd, dnegd 

720 write (i, 9915) denp, denn, denpd, dennd, atomp, atomn, atompd, 

& atomnd, epos, eneg, dif f (1) , dif f (2) , 

& temp, toler, zi (1) , zi (2) 

format (lx, 'MODEL INPUTS:'/ 


9913 
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0 

lx, ' ishow 

' , il5/ 

0 

lx, ' jmax 

' , il5/ 

0 

lx, ' nj 

' , il5/ 

0 

lx, ' nt 

' , il5/ 

0 

lx, ' ntpr 

' , il5) 

9914 format (lx, ' cur 

' , f 15 . 7/ 

0 

lx, ' area 

' , f 15 . 7/ 

0 

lx, ' cio 

' , f 15 . 7/ 

0 

lx, ' ecur (1) 

' , f 15 . 7/ 

0 

lx, ' ecur (2) 

' , f 15 . 7/ 

0 

lx, ' aa (1) 

' , f 15 . 7/ 

0 

lx, ' ac (1) 

' , f 15 . 7 / 

0 

lx, ' aa ( 2 ) 

' , f 15 . 7/ 

0 

lx, ' ac (2) 

' , f 15 . 7/ 

0 

lx, ' ocv 

' , f 15 . 7/ 

0 

lx, ' delta 

' , f 15 . 7 / 

0 

lx, ' dpos 

' , f 15 . 7/ 

0 

lx, ' dneg 

' , f 15 . 7/ 

0 

lx, ' dposd 

' , f 15 . 7 / 

0 

lx, ' dnegd 

' , f 15 . 7 ) 

9915 format (lx, ' denp 

' , f 15 . 7/ 

0 

lx, ' denn 

' , f 15 . 7/ 

0 

lx, ' denpd 

' , f 15 . 7 / 

0 

lx, ' dennd 

' , f 15 . 7/ 

0 

lx, ' atomp 

' , f 15 . 7/ 

0 

lx, ' atomn 

' , f 15 . 7/ 

0 

lx, ' atompd 

' , f 15 . 7 / 

0 

lx, ' atomnd 

' , f 15 . 7/ 

0 

lx, ' epos 

' , f 15 . 7/ 

0 

lx, ' eneg 

' , f 15 . 7/ 

0 

lx, ' dif f (1) 

' ,el5.7/ 

& 

lx, ' dif f (2) 

' ,el5.7/ 

& 

lx, ' temp 

' , f 15 . 7/ 

0 

lx, ' toler 

' , f 15 . 7/ 

0 

lx, ' zi ( 1 ) 

' , f 15 . 7/ 

0 

lx, ' zi (2) 

' , f 15 . 7/ 

0 

lx,' OUTPUTS:'/) 


c c(l,j) 

is the OH- concentration in mol/cm3 

c c (2, j) 

is the K+ concentration in mol/ cm3. 

c c (nspl. 

j) is the electric 

potential . 


ns=2 

nspl=ns+l 
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n=3 

del=delta/real (nj-1) 
two=2 . dO 
twodel=del*two 
f ar=964 87 . dO 

f rt=96487 . dO/ (8 . 314d0*temp) 
c 

C INITIAL CONDITIONS: 

DO 1 J=l, NJ 
c (1, j) =cio 
c (2, j ) =cio 
c (nspl, j) =0 .dO 
1 continue 

jcount=0 
jtd=0 
jtp=nt 
fmax=0 . dO 
time=0 . dO 
f lux=0 . dO 
f luxn j=0 . dO 


2 Calculate the initial amount of active material 

actpO=dpos*area*denp/atomp*epos*far 

actnO=dneg*area*denn/atomn*eneg*far 

actp(0) =actpO 
actn(O) =actnO 

tdel=dminl (dpos*0 . 999 d 0 *denp/atomp*epos, dneg*0 . 999d0* 

& denn/atomn*eneg) *far/cur ( jc) /real (nt) * (1 .dO-1 .d-6) 

aarea (1) =1 .dO 
aarea (2 ) =1 . dO 

c Print out initial conditions: 

write (2,293) j count, ocv, eta (1, 0) , eta (2, 0) , 

& aarea (1 ) , aarea (2) , cur ( jc) , time, actpO, actnO , 

& actp(0) /actpO, actn(0) /actnO 

write (2,2941) (real ( j-1) /real (nj-1) ,c(l,j),c(2,j), 

& j=l,nj) 

2941 format (3 (lx, ell . 5) ) 


c Start time loop: 
do 204 jt=l, nt 
write (6,722) jt,nt 

722 format (10X, 'Time step jt=',i3,' 
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c zk(jk) is used to time step the acid concentration, 

do 231 jk=l,nj 

231 zk ( jk) =c (1/ jk) 

time=real ( jt) *tdel 

C Calculate the amount of active material at t>0 
actp ( jt ) =actp ( jt-1 ) -cur ( jc) *area*tdel 
actn ( jt ) =actn ( j t — 1 ) -cur ( jc) *area*tdel 

aarea(l)=l .dO-time*cur ( jc) /dposd/denpd*atompd/epos/f ar 
aarea (2) =1 .dO-time*cur ( jc) / dnegd/ dennd*atomnd/eneg/f ar 
if (aarea ( 1 ). le . 0 ) go to 500 
if (aarea (2) . le . 0) go to 500 

c Find overpotential eta 
do 892 j=l, 2 
eta0=-2 . dl 
etal=2 . dl 

891 eta ( j, jt) = (etaO+etal) /two 

tcur=ecur (j) *aarea ( j ) * (dexp (eta ( j , jt) *f rt *aa ( j ) ) 

& dexp (-eta (j,jt)*frt*ac(j))) 

if (abs (tcur-cur ( jc) ) /cur ( jc) . le.toler) go to 892 
if (tcur . ge . cur ( jc) ) then 
etal=eta( j, jt) 
go to 891 
endif 

if (tcur . le .cur (jc) ) then 
eta0=eta ( j , jt ) 
goto 891 
endif 

892 continue 

c Start boundary value problem: 

24 JCOUNT= JCOUNT+1 
f lmax=0 . dO 
f 2max=0 . dO 
f 3max=0 .dO 

c Check number of program iterations: 

if ( j count .gt. jmax) then 
jtp= jt-1 

write (2,723) jt,nt 
write (3,723) jt,nt 
write (6,723) jt,nt 
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723 format (/, 10X, ' SOLUTION DOES NOT CONVERGE jt=' 

& ,13,' nt=',13) 

GOTO 500 

endif 

C Print to monitor current variable values: 
if (ishow . eq. 1) then 
write {*, 724) jcount , t ime, cur ( j c ) 

724 format (lx , ' jcount= ',115,' time 8 *' , e!5 . 7 , 

& ' cur=' , el5 . 7/ 

& ' zeta, cl, c2, c3,') 

write (*, 725) (real (js-1) /real (nj-1) ,c(l,js),c(2,js), 

& c (3, js) , js=l, n j) 

725 format (4 (lx, elO . 4 ) ) 

write (3, 7 241) jcount, time, cur ( jc) 

7241 format (lx, ' jcount= ',il5,' time=' , el5 . 7, 

& ' cur=' , el5 . 7/ 

& ' zeta, cl, c2, c3' ) 

write (3,7251) (real (js-1) /real (nj-1) ,c(l,js),c(2,js), 
& c (3, js) , js=l, n j) 

7251 format (4 (lx,el0 .4) ) 

endif 

do 2 j=l, n j 
do 2 k=l,n 

if ( jcount . eq • 1 ) z (k, j ) =c (k, j ) 

Z (k, j) = (1 .dO-dp) *c (k, j) +dp*z (k, j) 

2 continue 

c SET UP COEFFICIENT MATRIX 
J=0 

DO 31 1=1, N 
DO 31 K=l, N 
X (I, K) =0 . OdO 
31 Y ( I, K) =0 . OdO 
16 J=J+1 

DO 32 1=1, N 
DO 32 K=1,N 
G (I) =0 . OdO 
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A ( I, K) =0 . OdO 
B ( I, K) =0 . OdO 
32 D (I, K) =0 . OdO 

IF (J .NE. 1) GO TO 18 

C FIRST MESH POINT 

c Material-balance equation: 

fl=(-z(l,3)+4.d0*z (l,2)-3.d0*z(l,l))/twodel+ 
& cur ( jc) /zi (1) /far/diff (1) 

b (1, 1 ) =-3 .dO/twodel 
d(l,l)= 4 . dO/twodel 
x (1, 1) =-l .dO/twodel 

g(l)=-fl+b(l, 1) *z <1, j) +d(l, 1) *z (1, j+l)+ 

& x(l,l)*z(l, j+2) 

c Electroneutrality equation: 

f 2=0 . dO 

do 50 i=l,ns 

f 2=f 2+zi (i) *z (i, j) 

50 b (2, i) =zi (i) 

g(2)=-f2 

do 51 k=l,ns 

51 g(2)=g(2)+b(2,k)*z(k, j) 

c Potential equation: 

fnspl=z (nspl, j ) 
b (nspl, nspl) =1 . d0 

CALL BAND ( J) 

GO TO 16 

18 IF (J .EQ. NJ) GO TO 20 

C INNER MESH POINTS 

ck =zk(j) 
cl =z(l,j) 

del =(z(l, j+l)-z(l, j-1) )/twodel 

d2cl =(z (1, j+l)+z(l, j-1) -two*z (1, j) ) /del**2 

c2 =z (2 , j ) 

dc2 =(z(2, j+1) -z (2, j-1) ) /twodel 
c3 =z(3,j) 

dc3 =(z(3, j+1) -z (3, j-1) ) /twodel 

d2c3 = (z (3, j+1) +z (3, j-l)-two*z (3, j) ) /del**2 
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c Material-balance equation: 

f 1=- (cl-ck) /tdel+zi (1) *frt*diff (1) * 

& (cl*d2c3+dcl*dc3) +diff (1) *d2cl 

f icl=-l .dO /tdel+zi (1) *frt*dif f (1) *d2c3 
f lclp=zi (1) *frt*diff (1) *dc3 
f lclpp=dif f (1) 
flc3p=zi (1) *frt*dif f (1) *dcl 
f lc3pp=zi (1) *frt*diff (1) *cl 

c Band coefficients: 

a (1,1) =f lclp* (-1 • dO/twodel) + f lclpp* (1 .dO/del ) 

b ( 1, 1) =f lcl + f lclpp* (-two/del**2) 

d(l, i)=flclp* ( l.dO/twodel) + f lclpp* (1 • d0/del**2) 
a(l, 3)=flc3p* (-l.dO/twodel) + flc3pp* (1 .d0/del**2) 
b ( 1, 3) =f lc3pp* (-two/ del**2) 

d(l 3)=flc3p*( l.dO/twodel) + flc3pp* (1 .d0/del**2) 

gd) =-fl + a(l,l)*Z(l, j-l)+b(l,l)*Z(l, j) + d(l,l)* z (l, j + l) + 
& a(l,3)*z(3,j-l)+b(l,3)*z(3,j)+d(l,3)*z(3,j+l) 

c Electroneutrality equation: 

f 2=0 .dO 

do 52 i=l,ns 

f 2=f 2+zi ( i) *z (i, j) 

52 b(2 / i)=zi(i) 
g (2) =-f 2 

do 53 k=l,ns 

53 g(2)=g(2)+b(2,k)*z(k,j) 

c Potential equation: 

f 3 =dc3+ (cur ( jc) +far* (zi (2) *dif f (2) *dc2+ 

& zi (1) *diff (1) *dcl) ) / 

& (far*frt* (zi (2) **2*diff (2) *c2+ 

& zi (1) **2*dif f (1) *cl) ) 

f 3cl= - (cur ( jc) +far* (zi (2) *diff (2) *dc2+ 

& zi (1) *diff (1) *dcl) )/ 

& (far*frt* (zi (2) **2*diff (2) *c2 + 

& zi (1) **2*dif f (1) *cl) **2) * 

& zi (1) **2*dif f (1) 

f 3clp= zi (1) *diff (1) / 

& (frt* (zi (2) **2*dif f (2) *c2+zi (1) **2*diff (1) *cl) ) 

f 3c2= - (cur (jc)+far*(zi(2)*diff(2) *dc2+ 

& zi (1) *dif f (1) *dcl) )/ 
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& 

& 

& 

f 3c2p= 

& 


(f ar*f rt* (zi (2) **2*dif f (2) *c2 + 
zi (1) **2*dif f (1) *cl) **2) * 
zi (2) **2*dif f (2) 
zi (2) *dif f (2) / 

(frt*(zi(2)**2*diff (2) *c2+zi (1) 


**2*dif f (1) *cl) ) 


f 3c3p= l.dO 

a (3, l)=f3clp* (-l.dO/twodel) 
b (3, 1) =f 3cl 

d (3, 1) =f 3clp* ( l.dO/twodel) 
a (3, 2) =f 3c2p* (-l.dO/twodel) 
b (3, 2) =f3c2 

d(3, 2) =f 3c2p* ( l.dO/twodel) 
a (3, 3 )=f 3 c 3 pM -l.dO/twodel) 
d(3, 3) =f3c3p* ( l.dO/twodel) 
g ( 3 ) =-f 3+a (3, 1 ) *z (1, j-1) +b(3, 1) *z (1, 
; a(3, 2) *z (2, j-1) +b(3,2) *z (2, 

; a(3, 3) *z (3, j-1) 


j) +d (3, 1) *z (1, j+l) + 
j)+d(3,2)*z(2 / j+l) + 
+d (3, 3) *z (3, j + 1) 


Check on max f 1 , f 2, f 3, f 4, f 5 : 
f lmax=dmaxl (abs (f lmax ) , abs (f 1) ) 
f 2max=dmaxl (abs (f2max) , abs (f2) ) 
f 3max=dmaxl (abs (f3max) ,abs (f3) ) 


CALL BAND ( J) 
GO TO 16 


C LAST MESH POINT 


20 


del =UaJ" j-2)-4.d0*z(l, j-1) +3.d0*z (1, j))/twodel 

dc2 =U(2, 3 j-2)-4.dO*z(2, j-l)+3.d0*z(2 f j))/twodel 

c3 =z (3, j) _ 

dc3 = (z (3, j-2) -4 .d0*z (3, j-1) +3.d0*z (3, j) ) /twodel 


Material-balance equation: 

fl = dcl+ cur ( jc) /zi (1) /f ar/diff (1) 

y (1, 1) = l.dO/twodel 
a (1, 1 ) =-4 ,d0 /twodel 
b ( 1/ 1) = 3 . dO/twodel 

g ( 1 ) =-f 1+y (1,1) *z ( 1 , j-2) +a(l, l)*z (1, j-1) +b (1, 1) z (1, 3) 


c 


Electroneutrality equation: 
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f 2=0 . dO 

do 55 i=l,ns 

f 2=f 2+zi ( i) *z(i, j) 

55 b (2, i) =zi (i) 
g (2) =-f 2 

do 56 k=l,n 

56 g(2)=g(2)+b(2,k) *z(k, j) 


Potential equation: 

f 3 =dc3+ (cur ( jc) +far* (zi(2)*diff(2) *dc2+ 

& zi(l)*diff (l)*dcl) )/ 

& (far*frt* (zi (2) **2*diff (2) *c2+ 

& zi (1) **2*dif f (1) *cl) ) 

f 3cl= - (cur ( jc) +far* (zi ( 2 ) *dif f (2 ) *dc2+ 

& zi(l)*diff (l)*dcl) )/ 

& (far*frt* (zi (2) **2*diff (2) *c2 + 

& zi ( 1 ) **2*dif f (1) *cl) **2) * 

& zi (1) **2*dif f (1) 

f 3clp= zi (1) *diff (1) / 

& (frt* (zi (2) **2*diff (2) *c2+zi(l) **2*diff (1) *cl) ) 

f3c2= - (cur ( jc) +far* (zi (2) *dif f (2) *dc2+ 

& zi(l)*diff (l)*dcl) )/ 

& (far*frt* (zi (2) **2*dif f (2) *c2+ 

& zi (1) **2*dif f (1) *cl) **2) * 

& zi (2) **2*dif f (2) 

f 3c2p= zi (2) *diff (2) / 

& (frt* (zi (2) **2*diff (2) *c2 + zi (1) **2*dif f (1) *cl) ) 

f3c3p= l.dO 

y (3, 1) =f 3clp* (l.dO/twodel) 
a (3, l)=f3clp* (-two/del) 
b (3, I)=f3cl+f3clp*(3. dO/twodel) 
y (3, 2)=f3c2p* (l.dO/twodel) 
a (3, 2) =f 3c2p* (-two/del) 
b (3, 2) =f 3c2+f 3c2p* (3 . dO/twodel) 
y(3, 3) =f3c3p*( l.dO/twodel) 
a (3, 3) =f 3c3p* (-two/del) 
b (3, 3) =f 3c3p* (3 .dO/twodel) 

g (3) =-f 3 + y (3,l)*z(l,j-2)+a(3,l)*z(l / j-l)+b(3,l)*z(l,j) + 

& y(3,2) *z(2, j-2)+a(3,2) *z (2, j-1) +b(3,2) *z (2, j) + 
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y (3, 3) *z (3, j-2)+a(3,3)*z(3, j-1) +b (3, 3) *z (3, j) 


c Check on max fl,f2,f3: 

f lmax=dmaxl (abs (f lmax) , abs (f 1) ) 
f2max=dmaxl (abs (f2max) , abs (f2) ) 
f 3max=dmaxl (abs (f 3max) , abs ( f 3) ) 


CALL BAND ( J) 


C TEST FOR CONVERGENCE 


do 30 1=1,1 
do 30 j=l,nj 

if (z (i, j) .ne .0 .d0) change= dabs ( (c (i, j) -z (i, j) ) /z (i, j) ) 
if (c (i, j) . ne . 0 . dO ) change= dabs ( (c (i, j) -z (i, j) ) /c (i, j) ) 
if (c (i, j) .eq. z (i, j) ) change=0.d0 
if (change . ge . toler ) go to 24 

30 continue 

c List transient quantities: 

c Flux at j=l : 

c (The minimum of the coion or counterion flux determines 
c the acid or salt flux.) 

dcl= (-3 . d0*c ( 1 , 1 ) +4 .d0*c (1, 2) -c(l, 3) ) /twodel 
dc2= (-3 .d0*c (2, 1) +4.d0*c (2, 2) -c (2, 3) ) /twodel 
dc3= (-3 . d0*c (3,1) +4 .d0*c (3, 2) -c(3,3) ) /twodel 
f luxl=-zi (1) *frt*diff (1) *c(l, 1) *dc3-diff (l)*dcl 
f lux2=-zi (2) *frt*diff (2) *c (2, 1) *dc3-diff (2) *dc2 
flux=dminl (fluxl, flux2) 

c if (f luxl+f lux2 . le . 0 ) flux=dmaxl (fluxl, flux2 ) 

c Flux at j=nj: 

dcl=- (-3.d0*c (1, n j) +4 .d0*c (1, n j-1) -c (1, n j-2) ) /twodel 
dc2=- (~3.d0*c(2,nj)+4.d0*c(2, nj-1) -c(2, n j-2) ) /twodel 
dc3=- (- 3 .d 0 *c ( 3 , n j) +4 .d0*c (3, n j-1) -c (3, n j-2) ) /twodel 
fluxn jl=-zi (1) *frt*dif f (1) *c(l,nj) *dc3-diff (1) *dcl 
f luxn j2=-zi (2) *frt*diff (2) *c (2, nj) *dc3-diff (2) *dc2 
fluxn j=dminl (f luxnjl, f luxnj2) 
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if (fluxnjl+fluxn j2 . le . 0) f luxn j=dmaxl (f luxn jl, f luxn j2 ) 


fmax=dmaxl (abs (f lmax) , abs (f 2max) , abs (f3max) , abs (f 4max) , 
& abs ( f 5max ) ) 

c Potential at nj is the potential drop through membrane: 
pot ( jt) =c (3, n j ) 

volt ( jt)=ocv-eta (1, jt) -eta (2, jt) +pot ( jt) 
if (volt ( jt ) . ge . 0 . dO ) then 

if (actp ( jt ) . ge . 0 . dO ) then 

if (actn ( jt ) . ge . 0 . dO ) jtd=jt 

endif 

endif 

c Write results: 

if (amod (real ( jt) , real (ntpr) ) . It . 1 .d-9) then 
write (2,286) dcl,dc2,dc3, flux, f luxn j, pot ( jt) 

286 format ( lx,' ',/ 

0 lx,' del: ' , el2 . 5, ' 

0 ' dc3 : ' , el2 . 5/ 

0 lx,' flux: ',el2.5,' 

0 IX, ' IR DROP : ' , E12 . 5) 

write(2,293) jeount, volt ( jt) , 

& eta ( 1, jt) , eta (2, jt) ,aarea(l) ,aarea(2) , 

& cur ( jc) , time, actp ( jt) , 

& actn (jt) , actp ( jt) /actpO, actn ( jt) /actnO 

293 format ( lx, 'jeount ',il5/ 


& 

lx, ' volt 

' ,fl5.5/ 

& 

lx, ' eta (pos) 

' , f 15 . 5/ 

Sc 

lx, ' eta (neg) 

' , f 15 . 5/ 

Sc 

lx, ' aarea (pos) 

' , f 15 . 5/ 

Sc 

lx, ' aarea (neg) 

' , f 15 . 5/ 

Sc 

lx , 9 cur ( jc) 

' , f 15 . 5/ 

Sc 

lx , 9 time 

' , f 15 . 5/ 

Sc 

lx, ' actp ( jt ) 

' , f 15 . 5/ 

Sc 

lx, ' actn ( jt ) 

' , f 15 . 5/ 

Sc 

lx , 9 actp/actpO 

' , f 15 . 5/ 

Sc 

lx , 9 actn/actnO 

' ,fl5.5/ 

Sc 

lx,' 


Sc 

lx, ' j/n j, cl, c2 , c3' ) 



write (2, 294 ) (real ( j-1 ) /real (nj-l),c(l,j),c(2,j), 
& c (3, j) , j=l, n j) 

format (4 (lx, ell. 5)) 


dc2 : ' , el2 . 5, 

fluxn j : ' , el2 . 5/ 
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endif 

continue 


c Wirite results of transient quantities. 

500 power=volt (0) /two 
do 522 j=l, jtd 

522 power=power+volt ( j ) 

power=power*cur ( jc) /real ( jtd) / (delta+dpos+dneg) *1 . 0d3 
energy=power*t del* real ( jtd) /3 . 6d3 

cap =cur ( jc) *tdel* (real (jtd) +0 .5d0) / (delta+dpos+dneg) 
power l=power* (delta+dpos+dneg) 
energyl=powerl*tdel*real (jtd) /1000 .0d0 

capl =cur ( jc) *tdel* (real ( jtd) +0 . 5d0) 


2041 


526 


write (2,2041), cur ( jc) *1 . 0d3, cap, energy, power , 

& capl, energyl,powerl 

write (3, 2041) , cur ( jc) *1 . 0d3, cap, energy, power, 
i capl, energyl, powerl 

write (6, 2041) , cur ( jc) *1 . 0d3, cap, energy, power, 

£ capl, energyl, powerl 

format (/lx, 'Results of transient quantities:'/ 


lx, 'current density 
lx, ' capacity 
lx, 'energy density 
lx, 'power density 
lx, ' capacity 
lx, 'energy density 
lx, 'power density 


, f 15 . 5, ' 
, f 15 . 5, ' 
' ,fl5.5,' 
' , f 15 . 5, ' 
' ,fl5.5, ' 
' ,fl5.5, ' 
' , f 15 . 5, ' 


mA/cm2 ' / 
kC/1' / 
Wh/1' / 

W /l'/ 

C /cm2' / 
Ws/cm2' / 
mW/cm2 ' / 


lx, ' capacity, potential') 


do 527 j=2, 3 

write ( j , 2042 ) 0.d0,ocv 

write (j, 2042) (cap*real ( jt) /real ( jtd) , 

volt ( jt) , jt=ntpr, jtp,ntpr) 


527 write (j, 2043) 

write (6, 2042) 0.d0,ocv 

write (6, 2042) (cap*real ( jt) /real (jtd) , 

& volt ( jt) , jt=ntpr, jtp, ntpr) 

2042 format (2 (lx, elO . 4) ) 

2043 format (//) 

501 CONTINUE 

502 CONTINUE 

503 CONTINUE 

504 CONTINUE 

505 CONTINUE 
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stop 

end 


SUBROUTINE MATINV (N, M, DETERM) 
implicit double precision (a-h,o-z) 
COMMON A(7, 7) ,B (7, 7) , C (7, 500) , D (7, 15) 
DIMENSION ID (7) 

DETERM=1 . 0 
DO 1 1=1, N 
1 ID (I) =0 

DO 18 NN=1 , N 
BMAX=1 . 1 
DO 6 1=1, N 

IF {ID (I) . NE . 0 ) GOTO 6 
BNEXT=0 . 0 
BTRY=0 . 0 
DO 5 J=l, N 

IF (ID ( J) . NE . 0 ) GOTO 5 

IF (dabs (B (I, J) ) .LE.BNEXT) GOTO 5 

BNEXT=dabs (B (I, J) ) 

IF (BNEXT.LE.BTRY) GOTO 5 

BNEXT=BTRY 

BTRY=dabs (B (I, J) ) 

JC=J 

5 CONTINUE 

IF (BNEXT . GE . BMAX*BTRY) GOTO 6 

BMAX=BNEXT/BTRY 

IROW=I 

JCOL= JC 

6 CONTINUE 

IF(ID(JC) . EQ . 0 ) GOTO 8 

DETERM=0 . 0 

RETURN 

8 ID ( JCOL) =1 

IF ( JCOL . EQ . IROW) GOTO 12 
DO 10 J=1 , N 
SAVE=B (IROW, J) 

B ( IROW, J) =B { JCOL, J) 

10 B ( JCOL, J) =SAVE 
DO 11 K=1 , M 
SAVE=D (IROW, K) 

D (IROW, K) =D (JCOL, K) 
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11 D ( JCOL, K) =SAVE 

12 F=1 .0D0/B( JCOL, JCOL) 

DO 13 J=1 , N 

13 B (JCOL, J) =B (JCOL, J) *F 
DO 14 K=1,M 

14 D (JCOL, K)=D (JCOL, K) *F 
DO 18 1=1, N 

IF (I. EQ. JCOL) GO TO 18 
F=B(I, JCOL) 

DO 16 J=1 , N 

16 B ( I, J) =B (I, J) -F*B (JCOL, J) 

DO 17 K=1 , M 

17 D (I, K) =D (I, K) -F*D (JCOL, K) 

18 CONTINUE 
RETURN 
END 

C 

C 

C 

SUBROUTINE BAND ( J) 

implicit double precision (a-h,o-z) 
DIMENSION E (7, 8, 500) 

COMMON A (7, 7) ,B(7, 7) ,C(7, 500) , D (7, 15) , 
& G(7),X(7,7),Y(7,7),N,NJ 

101 FORMAT (15H DETERM=0 AT J=,I4) 

IF ( J-2) 1,6,8 

1 NP1 = N + 1 
DO 2 1=1, N 

D ( I, 2*N+1 ) = G (I) 

DO 2 L=l, N 
LPN = L + N 

2 D (I, LPN) = X (I, L) 

CALL MATINV(N, 2*N+1, DETERM) 

IF (DETERM) 4,3,4 

3 write (3,101) J 

4 DO 5 K=l, N 

E (K, NP1, 1 ) = D (K, 2*N+1 ) 

DO 5 L=l, N 
E (K, L, 1 ) = - D (K, L) 

LPN = L + N 

5 X (K, L) = - D (K, LPN) 

RETURN 

6 DO 7 1=1, N 
DO 7 K=l, N 
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DO 7 L=1,N 

7 D ( I, K) = D ( I , K) + A ( I , L) *X (L, K) 

8 IF (J-NJ) 11,9,9 

9 DO 10 1=1, N 
DO 10 L=1,N 

G ( I) = G ( I ) - Y(I,L)*E(L,NP1, J-2 ) 

DO 10 M=1 , N 

10 A ( I, L) = A (I, L) + Y(I,M)*E(M,L, J-2) 

11 DO 12 1=1, N 

D ( I, NP1 ) = - G ( I ) 

DO 12 L=l, N 

D ( I, NP1 ) = D ( I , NP1 ) + A(I,L) *E(L,NP1, J-l) 

DO 12 K=1 , N 

12 B (I, K) = B ( I , K) + A (I , L) *E (L,K, J-l) 

CALL MATINV(N,NP1, DETERM) 

IF (DETERM) 14,13,14 

13 write (3,101) J 

14 DO 15 K=1 , N 
DO 15 M=l, NP1 

15 E (K, M, J) = - D (K, M) 

IF (J-NJ) 20,16,16 

16 DO 17 K=1 , N 

17 C (K, J) = E (K, NP1, J) 

DO 18 JJ=2, NJ 

M = NJ - JJ + 1 
DO 18 K=1 , N 
C (K, M) = E (K, NP 1 , M) 

DO 18 L=l, N 

18 C (K, M) = C (K, M) + E(K,L,M) *C(L, M+l ) 

DO 19 L=1 , N 

DO 19 K=1 , N 

19 C (K, 1 ) = C (K, 1 ) + X (K, L) *C (L, 3) 

20 RETURN 
END 
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